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SUMMARY 


The present analysis is developed to calculate the com- 
pressible three dimensional inviscid flow in the scroll and 
vaneless nozzle of radial inflow turbines. A Fortran computer 
program which was developed for the numerical solution of this 
complex flow field using the finite element method is presented. 

The program input consists of the mass flow rate and 
stagnation conditions at the scroll inlet and of the finite 
element discretization parameters and nodal coordinates. The 
output includes the pressure, Mach number and velocity magnitude 
and direction at all the nodal points. 

The report includes the computer program, the analysis, 
a description of the solution procedure, and a numerical 
example . 
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INTRODUCTION 


Existing quasi-three-dimensional turboraachinery flow field 
solutions are applicable to radial inflow turbine rotors but 
not to its scroll. Most analytical studies of the flow in the 
scroll [1, 2] are based on the assimption of one dimensional 
flow and the conservation of mass and angular momentum. A 
two dimensional study of the flow field in the scroll and 
nozzles was reported in reference [3] . Another two dimensional 
study of the velocity components in the cross-sectional planes 
of the scroll was reported in reference [4] . These two 
analytical investigations provided insight into the influence 
of the scroll and nozzle vane geometries on the flow field 
but together do not constitute a quasi-three-dimensional 
solution. Experimental measurements of the velocity components 
in the radial turbine scroll [5] indicate that actually 
the three dimensional scroll flow fields are very complex. 

This work presents a method for the analysis and a 
code for computing the compressible inviscid three dimensional 
flow field in the scroll and vaneless nozzle of a radial 
inflow turbine. The finite element method is used in the 
numerical solution in order to model the complex geometrical 
shape of the flow boundaries. 

MATHEMATICAL ANALYSIS 

The three basic steps involved in the finite element method 
are the domain disretization, obtaining the integral statement 
of the problem and the appropriate choice of the local approxi- 
mation in each element [6, 7]. In this study, the three 
dimensional flow field is formulated in terms of the potential 
function and Galerkin's method is used to obtain the integral 
statement. 

For steady flow, the equation of conservation of mass is 
expressed in terms of the potential function (}> as follows; 
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V* (pV(})) = 0 


( 1 ) 


The unknown velocity potential solution is approximated in 
each element by 

4)® = [N] (2) 

where is the column vector of the potential function 

values at the nodes of the element e, and [N] is the row 

vector of the interpolation functions 

Galerkin's method requires the residual of the substitution 
of equation (2) into equation (1) to be orthogonal to the 
interpolation functions 

/ N. V (pV4)) dV = 0 . (3) 

V ^ 

Applying Green's theorem to equation (3) leads to the following 
relation which involves the boundary conditions : 

/ VN. -(pVtJ)) dV = / N. pV(f)*dA (4) 

V ^ A ^ 

The integrals in the above equations are evaluated over the 
volume V and surface A of all the solution domain. The right- 
hand side of equation (4) represents the weighted average of 
the mass flow across the solution domain boundaries. The solid 
boundaries contribution to equation (3) is equal to zero because 
of the no-flux condition. 

Equation (4) holds for the entire solution domain and also 
for an arbitrary finite element e. The subsntution of equation 
(2) into equation (4) gives the following element equation: 

/ VN^-(pV[N] {())}®)dV = / p dA (5) 

V® A® 

Equation (5) is nonlinear for compressible flow where the 
density is given by; 


3 



( 6 ) 


P 



RT ^ 
o 


-1/Y-l 


The numerical solution to equations (5) and (6) is obtained 
iteratively by allowing the density, p, in equation (5) to lag 
one iteration behind the potential flow field solution. The 
assembly of all the linearized element equations in the flow 
field results in a set of linear algebraic equations which can 
be expressed as follows; 

[K] {4>} = {P} (7) 

The Gauss elimination method is used in the solution of the 
linearized equations. 

The matrix of coefficients is updated at each iteration 
but the load vector remains unchanged, as it is determined by 
the mass flow rate in the scroll. Appendix A gives a detailed 
derivation of the element equations and of the coefficients 
of the stiffness matrix. 


The Splitting Surfaces 

The choice of the potential function as a dependent variable 
is generally suitable for three dimensional flow field solutions 
in simply connected domains. However the domain in the scroll 
problem is multi-connected leading to a mathematically non- 
unique solution for the potential function. In the present 
analysis, two coincident surfaces are introduced in the solution 
domain to make it simply connected. The splitting surfaces 
extend from the scroll lip, which is the end of the metal between 
the inlet section and the last scroll section, to the vaneless 
nozzle exit station as shown in Fig. 1. In the flow field 
discretization, two sets of nodes are placed on the opposite 
sides of the splitting surfaces and a discontinuity in the velocity 
potential field is allowed across. Additional matching conditions 
are therefore required to insure the continuity in the velocity 
field across these boundaries as described in Appendix B. 
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INPUT 


The input data to the program is read in subroutine INPUT. 

It consists of the geometrical data and the performance para- 
meters for a given scroll. Any consistent set of units can be 
used for the input variables. Other input parameters are related 
to the numerical solution procedures and to produce the desired 
output. Following is a description of the input data in 
alphabetical order. 


BANDWD 

GAMA 

EPSILN 

INDATA 

ITERMX 

J1(I) 


An integer equal to the band width of the matrix 
of coefficients. 

Specific heat ratio. 

A small number describing the tolerance which is 
used to determine the convergence of the numerical 
solution based on the normalized density change. 

Integer to indicate if it is desired to print out 
the input data. INDATA = 0 will not produce a 
printout of the input data while INDATA = 1 will 
produce a printout of the input data. 

The number of iterations not to be exceeded in 
the solution. Incompressible flow solution is 
obtained if ITERMX = 1. 

An array of the global number of the first node 
of every cross section I 


LBC(I) ,MBC(I) , 
NBC (I) 

MS RATE 

MSI 

MS2 

MS3 


Arrays of the nodal points of the surface element I. 

The mass flow rate in the scroll (kg/sec or slg/sec) . 

The nvimber of surface elements at the inlet station. 

Total number of surface elements at inlet section 
and the first splitting surface. 

Total number of surface elements at inlet section 
and the two splitting boundaries. 
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MSTOT 

MTOT 

MV(I) 

NANZ 

NC 

NRNZ 

NODEd, J) 

NTOT 

NU(I) 

R 

STAG? 
STAGT 
THETA (I) 


X(I) ,Y(I) , 
Z(I) 


Total number of surface elements at inlet, exit 
and splitting boundaries. 

Total number of volume elements . 

The global number of the volume element which con- 
tains the splitting surface element I (I = MSl+1 to MSS). 

The number of nodal points in the axial direction 
in the nozzle. 

Total number of cross sections in the scroll and 
inlet passage. 

The number of nodal points in the radial direction 
in the nozzle. 

Array of the global number of the nodes, 

I = 1,2, 3, 4, of the element J. J = 1, MTOT. 

Total number of nodes. 

An array of the number of nodes at every cross 
section I. 

Gas constant: (J/°K or lb ft/slg °R) . 

2 3 

Stagnation pressure (newton/m or Ib/ft ) . 

Stagnation temperature (°K or °R) . 

An array of the angles between every cross section 
plane and the x-z plane measured from the positive 
x-direction. 

Arrays of the Cartesian coordinates of all the 
nodal points. 


The reading foinnats of the input variables are illustrated 
in Table 1. 


Input Preparation 

The solution region includes the inlet, the scroll and the 
vaneless nozzle. The splitting boundaries should extend from 
the lip (the last metal part between the inlet and the last 
scroll section) to the exit station of the vaneless nozzles. 
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It can have any shape that is consistent with the discretization 
pattern. The discretization pattern and system topology is 
arbitrary. The surface elements at the inlet, exit and the 
splitting boundaries are inputted in one array. They are 
numbered sequentially up to MSI on the inlet section, then up 
to MS 2, then MS3 on the two splitting boundary surfaces, 
then up to MSTOT on the exit station. The volume elements 
associated with these surface elements are inputted in the 
array MV (I), 1=1, MSTOT. The angular coordinates THET(I), 
of the cross sectional planes, including those at inlet, are 
only used to evaluate the through flow velocity and the 
radial velocity from the computed velocity components in 
Cartesian coordinates after the convergence of the numerical 
solution. 

It is up to the user to select the number of cross sections 
and the nodal points per cross section. The band width is 
equal to the maximum difference between the global number of 
any two nodes of the same element. The system topology can be 
specified in such a manner as to minimize the band width of 
the matrix of coefficients in the numerical solution. For the 
purpose of organizing the printout of the output, the nodes 
in each section N(J) are numbered sequentially starting with 
JI (J) , J = 1,NC. 


OUTPUT 

A sample output is given in Table 2 for the scroll of 
the numerical example. The first set of output consists of a 
printout of the input parameters. The second set of output 
consists of a listing of the potential function, the through 
flow velocity, the cross velocities, the total velocity, the 
Mach number, the pressure and density normalized by their inlet 
stagnation values. The computed flow data are listed for 
every section starting at inlet so that they can be easily 
interpreted. 
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Error Messages 

The following error conditions lead to program termination 

1. ELEMENT NUMBER XXX HAS A ZERO VOLUME 
(LESS OR EQUAL TO XXX) 

CALCULATION IS STOPPED. 

The message is printed in subroutine MATVEC during the first 
iteration as the element volumes are evaluated in ascending 
element order. Check the coordinates of NODE (I; ELEMENT 
NUMBER XXX) for error. 

2. WARNING: 

SINGULARITY IN ROW XXX 
CALCULATION IS STOPPED. 

This message is printed in subroutine GAUSS if the diagonal 
element of the matrix of coefficient is less than 10 The 

row number is printed to be used for the identification of the 
corresponding nodal point. Check the topology of the elements 
associated with this node in the discretization pattern. 

PROGRAM DESCRIPTION 

The main program is designed to obtain the potential flow 
solution in a radial inflow turbine scroll using the finite 
element technique. The program is built from four major 
subroutines: MATVEC, SPLIT, GAUSS and PHYSIC, two edit sub- 
routines: INPUT and OUTPUT, and four supplementary subroutines 

VOLELM, NUMBER, AR and COFACT which are called inside the 
major subroutines at various steps for side calculations. 

The structure of the program is shown on Fig. 2. 

Descriptions of the Subroutines 
I. Major subroutines: 

MATVEC: The stiffness matrix elements, , and the load 

vector, P j , are calculated in this subroutine. In addition, 
the volume of the elements are calculated and stored in the 
array VOL (I) so that it can be used in other subroutines and 
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in subsequent iterations. The coefficients of the element stiff- 
ness matrix, are calculated using Eq. (A. 8). Subroutine COFACT 

is called to calculate the coefficients a, , b, and c, (Eq. A. 7). 

K jC K 

The load vector components are calculated using Eqs. (A. 10). 

SPLIT: The matching conditions across the splitting surfaces, 

as explained in Appendix B are implemented in this subroutine. 

GAUSS : The solution for the potential function at the 

nodal points (I) is obtained in this subroutine using the GauSs 
elimination method. The matrix of coefficient is bounded, but 
nonsymraetric due to the application of the matching condition at 
the splitting boundary. Since the solution to the potential 
function is unique within an arbitrary constant, the magnitude 
of the potential function at the first node is arbitrarily 
fixed at -1. 

PHYSIC: The physical properties of the fluid are calculated 

in this subroutine which is called after the nodal potentials 
are determined. The properties calculated at each node are: 
the velocity components as defined by Eq. (A. ) , the Mach 

number and the normalized pressure and density. The difference 
between the element density of two subsequent iterations are 
also evaluated in this subroutine to be used in the convergence 
criteria. 

II. Edit subroutines: 

INPUT: Called first in SOLVE to read and write the input 

data. It also initializes the normalized density in each 
element to unity. 

OUTPUT: After the convergence of the solution, the velocity 

potential and flow properties are printed at all the nodes. The 
printout is organized so that the data of each cross section 
is identified and printed separately. 
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III. side calculations subroutines; 

NUMBER is called by SPLIT to rearrange, at the element level, 
the nodes of the elements with surfaces on splitting boundary, 
such that the first three nodes of an element are on the splitting 
boundary. The nodes global numbers are then stored in the 
array NSB(I,J), where I is the number of the surface element 
on the splitting boundary and J = 1,2, 3, 4 are the global numbers 
of the three nodes on the splitting boundaries followed by the 
fourth internal node.. 

VOLELM is a subroutine with double indices N and K. The 
first index N refers to the element, and the second is a flag. 

If the index N is equal to 2, the coefficients of the interpolation 
functions are calculated according to equation (A. 7). If it is 
equal to 1, a similar procedure is followed for the voliame 
elements which have three nodes on the splitting surface whose 
nodes are locally renumbered and identified by the number of the 
surface element. If the second index is equal to zero, the 
volume of the element is also calculated according to equation (A.9) . 
The volume calculations are only performed during the first 
iteration, then stored in VOL (I) for subsequent use in later 
iterations. 

IV. Dictionary of common block variables; 

A(I,J) Global matrix of coefficients for the potential 

function at the nodes. 

ART (I) Area of surface element I (at inlet, exit and/or 

splitting boundary) . 

B(I,J) The matrix of the element nodal coordinates. 

E(I) The contribution of mass flow rate through a splitting 

boundary surface element, to the opposite surface 
element equation, in the global matrix of coefficients. 

JI(J) The global number for the first node in a scroll 

cross section J. 


10 



LBC(I) ,MBC(I) , 
NBC ( I ) 

LOAD (I) 

MV(I) 

NU(J) 

NODE(I,K) 

NTL ( I ) 

NSB(I,J) 
PRTIALd, J) 

VOL ( I ) 

X(I) , Y(I) , 
Z(I) 


The nodal number for the three nodes of a 
surface element. 

An array containing the load vector. 

The global number of a volume element which con- 
tains the splitting boundary surface element I. 

Total number of nodes in the scroll cross 
section J. 

The global number for the node of the voliame 
element K, for the tetrahedral volume element 
I = 1,2, 3, 4. 

An array containing the nodes on a splitting 
boundary arranged in descending order. 

The global number of the nodes I of the element J 
with 3 nodes on the splitting surface. 

Partial derivatives of the interpolation functions 
Nj (J = 1,2, 3, 4) in Cartesian coordinates, 

IsF ^ ly I = 2, for I = 3) . 

Volume of an element I. 

Coordinates of a node I. 
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PROGRAM LISTING 


12 



nnnnnnononon n n nnonnn non 


AN IV Gl RELEASE 2.0 


MA IN 


DA TE 


82239 ' 


04/25/19 


INTEGER BANDWD 

REAL MSRATE .MSFLUX .L0A0.MACH2 

COMM ON/ASEMOL/H( 4,4) . VOL ( 3720 ), NODE ( 4 , 3720 ) ,A(9t> 0,406) ,NTL(155) , 

. PHTI AL(3,4), ART(155),LHC( 155),MRC( 1 55 ) ,N BC ( 1 55 ) , 

. X(960) ,y<960) ,Z(960) , LOAD (960) ,nANOWO,E( I 55), MV ( 155) 

. , NS8( 4, 1 55) , MSI , MS2 , MSS , MSTOT , NTOT , JUMP 1 , JUMP2 ,NANZ , 

. NRN2,POTEN(960 ) 

COMMON/PHYSI /G AM A, R , ST AG T , ERROR , VELO ( 3 ) 

COMMUN/INOUT/INDATA,MTOT,MSRATE,EPSILN,I TER, I TERMX , NC , ST AGRO , 

, Jl (60 ) ,NU(60 ) , THETA ( 60 ) ,ROELEM( 3720 ) , NN 

» WARNING; THE DIMENSION OF A(I,J) AS IN COMMON SHOULD AT LEAST 

BE J=( RANDWO-1 ) *2 . 


CALL UNDFLW 

CALL INPUT 
(TER = 0 : 

NN = BANDWD«2 - 1 
DO 80 I=l,NTOT 
aOl LOAD! I ) = 0.0 

10 ITER = ITER + I 
DO 90 I=l,NTOT 
DO 90 J=1 ,NN 
90 A( I , J) - 0.0 

BUILD UP THE STIFNESS MATRIX AND THE LOAD VECTUP. 
CALL MATVEC 


INTRODUCE THE SPLITTING SUFACE MATCHll'lG CONDI FIUNS. 
CALL SPLIT 


SFT THE rinST DIAGONAL TERM IN STIFNESS flATRIX AND ITS 
CUKRESPONDINNG VALUE IN LOAD VECTOR EQUAL TO UNITY. 

K - OANDWD + 1 

r = ABS( A ( 1 .OANOWO ) ) 

DO 160 I=K,NN 



on non nnn nnn nnnn non 


IV Gl RELEASE 2.0 


MA IN 


DATE 


Q22J0 


04/25/19 


160 A(l*n = I .0E-15*A(l )/F 

A(1,0ANDWO) = A{ UBANOWD )/F 
lOAD(l) = LOADt I )/ABS(L0AD( 1 > ) 

IflPOSE ZERO TANGENTIAL VELOCITY AT INLET SECTION. 

NODIN = NU( 1) 

DO 170 I =2, NODiri 
L = K - I 

A( I .L) = A( I .L) + 1 . 0E4^15 
170 A(I.BANl)WD) ^ A(I*BANDW0) - l.OE-MS 

SEAK A SOLUTION FOR VELOCITY POTENTIALS; THE INLET SECTION IS 
AN EOUIPOTLNTIAL OF MAGNITUDE UNITY. 

CALL GAUSS 

CALCULATE THE PHYSICAL PROPERTIES OF THE FLUID. 

CALL PHYSIC 

tEST THE CONVERGENCE AND ITERATE. 

IF (ERROR. GT.EPSILN .AND. ITER .LT . ITERMX ) GO TO 10 

PRINTOUT THE CALCULATED PHYSICAL ROPERTIES IN THE CURRENT 
IT£F^ATI3N AT ALL NODES: VELOCl TY PUTENTIALS. VELOCITY COMPONENTS# 

MACH NUMBER NORMALIZED PRESSURE AND DENSITY. 

CALL OUTPUT 

5>T0P 

END 



\N IV G1 RfLFASC 2.0 


MATVEC 


DATE = 82239 04/25/19 


Ul 


C 

c 

c 

c 


c 

c 

c 


c 

c 

c 


c. 


c 

c 

c 

c 


c 

c 

c 

c 


SUOnOUTINE , MATVEC 

THIS SUDROUTINE BiUlLOS UP THE STIFNESS MATRIX AUO THE LOAD 
VECTOR * 

TraTEGER BANDWO 

HEAL MSRATE.MSFUUXtLOADtMACH^ 

COMMON/ ASEMBL/BC 4,4) . VOL ( 37 20 ) • NODE ( 4* 3720) » A( 960,4 0 6) #NTL( 155) . 

. PRTI AL(3 .4 ) .ART (155 ) .LBC( 155 ) .MF)C( 1 55 ) , N BC ( 1 55 ) , 

. X( 96 0) , Y{960) .7(960) .LOAD( 960 > . BANDWO. E( 155 ) .MV( 155 ) 

. . NS 8 (4» 155) .MS1.MS2.MS3.MSTOT.NTOT. JUMPl. JUMP2.NANZ, 

. NRN2,POTENC960) 

COMMON/INQUT/INDATA.MTOT.MSR ATE. EPS! LN, I TER . I TER MX . NC . ST AGRO , 

. J1 (60 ) .NU (60 ) » THETA ( 60 ) . ROELEMC 37 20 ) . NN 


I FORMAT ( IHl .20(/) .25X. 'ELEMENT NUMBER' . I 1 5 . / . 2 5 X . 

. 'HAS A ZERO VOLUME (LESS OR EQUAL Ta*.lPF|2.4. 

. •)'. /.25X. • CALCULATION IS STORED.') 


MARY = 0 

IF (ITER .NE. 1) MARY = 2 

DO 120 J[ 1= 1 » MTOT 
CALL VOLELM(MARY. I I ) 

IF (VaL(II) .GT. l.OE-12) GO TO 100 
II .VOL( I I ) 

STOP 


CALCULATE THE COEFICIENTS OF X . Y, & Z IN SHAPE FUNCTIONS AND THE ‘ 
PARTIAL DERIVATIVES OF THE SHAPE FUNCTIONS. 


I 00 


DO 110 1=1 .4 

I 10 J= I . 3 
Pk r 1 AL ( J . I ) = 


COFACTI I . .B ) 


CALCULATE THE ELFMENT 
GLOBAL STIFNESS MATRIX 


STIFNESS MATRIX 
A( K ,i ) . 


•ELEM' AND BUILD UP THE 


F = 

i ')0 

DU 


ROELEM( I I )/( VOL( I I ) ^ 3b ) 
120 1=1,4 

I 20 J = 1 .4 


L LEM 


( PRTI AL ( 1 . I )*PRTI AL( I , J ) + 
Pk T I AL( 2, I ) <^PR r I AL ( 2, J ) + 
PRTI AL( 3 , I ) <^PPT I AL( 3 , J) ) 


I 10 



AN IV G1 RELKASE 2.0 


MATVEC 


DATE = Q2->39 


04/25/19 


K = NOOEC I . I I) 

L = NODE(J.II) - NODE! I, II) *■ BAt^DWO 
120 A{K.L) = A(K,L) + ELEM 

IF (ITER .GT. 1) GO TO 170 

C 

C BUILD UP THE LOAD VECTOR LOAD(I). 

C 

, AREAIN = 0.0 

AHEAXT = 0.0 
C 

DO 130 I=1>MS1 

' r ART( I ) = AR( I ) 

, 130 AREAIN = AREAIN + ART(I) 

C 

MS =: MS3 + 1 
DO 140 r=MS«MSTOT 
ART(I) = AR(I) 

140 AREAXT = AREAXT ART(I) 

C 

M MSFLUX = MSRATE/I AREAIN*STAGRO) 

DO 150 I^l.MSl 
P = - MSFLUX+ARTI I 1/3 
LOAD(LBC(in = LOADILBCID) + P 

LOADIMBCd)) = LOAD(MBC(I)) + P 

ISO LOADINBCII)) = LOAD(NPC(I)) + P 

C 

MSFLUX = MSRATE/( AREAXT»STAGRO) 

DO 160 I=MS,MSTOT 
P = MSFLUX* ART ( I 1/3 
LOAIXLRC ( I ) ) = LOAD(Lf!C( I ) ) 4 - P 

LOAD (MBC ( I ) ) = LOAD(MDCd)) + P 

160 LOADINBC(I)) = LOADINBC(I)) F P 

C 

< 170 DO 180 I=l,NTOT / 

100 POTtN(I) = LQADII) 


C 


RETURN 

END 
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AN IV G1 RELFASH 2.0 


SPLIT 


DATF- 


8 2239 


04/25/19 


SUBROUTINE SPLIT 

THIS SUOROUTINE MODIFIES THE STIFNESS MATRIX BY INCLUDING THE 
SPLITING BOUNDARY CONTRIBUTION. 


INTEGER BANDWD 

REAL MSRATE, MSFLUX .L0A0.MACH2 

CaMM0N/ASEM8L/D(4 .4 ) .VOL I 3720) * NODE I 4, 3720 ) .A(96 0,4 06) ,NTL( 155 ) , 

. PRTIALI 3. 4). ART! 155) ,LBC( I 55 ) , MRC( 155) .NBC (155) , 

. X ( 960 ), Y ( 960 ). Z (960), LOAD! 960), BANDWD, E( 155) ,MV( 155) 

. , NSB( 4, 1 55) ,MS1 ,MS2 , MS3 , MSTOT ,NTOT , JUMP 1 , JUMP2.NANZ, 

, NRNZ,POTEN(960) 

COMMON/INOUT/INDATA,MTOT, MSRATE, EPS ILN. ITER, I T ERMX , NC . ST AGRO . 

. Jl(60),NU(60) , THETA ( 60 ) ,ROELEM(1720 ) , NN 


IF (ITER .EQ. 1) CALL NUMOER , 

START THE MAIN LOOP OVER THE VOLUME (OR, SURFACF) ELEMENT. 

/ 

MDir = MSI - MS2 
MS = MSI + 1 
MSS = MS2 
KLM = MS 

‘ \ 

10 DO 60 II =MS,MSS 

BUll/D UP THE MATRIX B(I,J). 

CALL VOLELM( 1,11) 

CALCULATE THE PARTIAL DERIVATIVES Of SHAPE FUNCTIONS ' 

DO 20 1=1.4 

00 20 J=1.3 

20 PRTIAL(J.I) = COFACT( I , J+l ,B ) 

CALCULATE THE AREA OF TfiE SURFACE ELEMENT II. 

ARFA = SORT (PRTI AL( 1 .4 )**2 + PRT I AL ( ? , 4 ) +♦ 2' I PR T I AL ( 3 , 4 ) ♦ ♦ 2 ) / 2 
IF (II.NE.KLM .OR. ITER.NE.l) GO TO 30 

EVALUATE TEIE COEFFICIENTS OF THE PLANE EQUATION A*X E B+Y + 

Cifl + D = 0 AND THE DIRECTIONAL COSINES OF THE COMPONENTS OF 
THE NORMAL VECTOR. 

AX = PRTIALI 1,4) 
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BY = PRTIALC2.4) 

, ■ , ' ^ CZ ^ PRTI ALI3.4) ' . , 

' . . O =r B( 1*2)4(B(3« 3)«B(2*4) - B( 2* 3) 4Q( 3 «4 ) ) 4- 

' • B(l.3)*(B(3t4J*B(2,2) - B<2.4)*B(3.2 ) ) 4^ 

_ • B(li4) 4(B(3»2)«B(2.3^ - B( 2*2)4B(3.3) ) 

DELTA - AREA+2 ’ ■ * / ’ ' < ' , ' > 

->■ ■ IF (D. ME. 0*0 *AND. O^DELTA^GT. 0. 0 ) DELTA = -DELTA 

, ' IF (D«EQ*0*0 *AND, CZ.NE.O.O .AND. CZ*OELTA.LT. 0. 0» DELTA = -DELTA 

■ IF IO.EQ.0.0 .AND., CZ.EQ.O .0 .AND. BY+DELTA ,LT .0 .0 I ' DELTA = -DELTA 
COSINA = AX/DELTA , 

COSINB = DY/OELTA - < ' . > 

CQSINC = CZ/DELfA. ^ 

, EVALUATE THE COEFFICIENTS Et I » OF fTHE NODAL POTENTIALS IN A | 

BOUNDARY ELEMENT ,I I . ' ' 

T - t * , ^ ' 

30, F= ROELEMIMVIIDI/CVOLIMVI II))*18) ’ , ■ '' 

DO 40 l=l .4 

Eli) = AREA»tCOSINA*PRTI ALI 1 .1) + C0SINB4PRTI AL ( 2. I) + 

. COSINC4PRTIALI3. 1 ),)4F 

40 IF (MSS .EQ.: MS3) Ed) - -E ( I ) 

/ ^ 

INSERT THE SURFACE INTEGRAL IN STIFNESS MATRIX. * 

IJ ^ II - MOIF ; 

, DO 50 J=l , 4 
DO 50 I =1 ,3 
K = NSD( 1,11) 

L = NSB(J,II) - NSBd.IJ) + BANDWD 
50 A(K.L) = AIK.L) + E(J) 

60 CONT INUE 

MDIF = -MDIF ' 1 

MS = MS2 + I ^ 

MSS = MS3 , ,'i 

IF (MDIF .GT.-O) GO TO 10 

INCLUDE TENGENTIAL VELOCITY IN STIFNESS MATRIX TO INSURE 
CONTINUITY THROUGH SPLITINQ BOUNDARIES. t 

JUMP = JUMP! ' I 

K1 = NTL{ I ) 

K2 = K1 - JUMP 
MSO = (MST - MS2)»3 
MONA = Kl + I -( NRNZ- 1 ) 4NAN2 
C 

DO 80 I=2,MSD 
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IF (NTLtl) .EQ. NTLM-Dl GQ TO 80 
i IF t NTLtl) .GE. MONA) GO TO 70 

RETURN 

I 65 K2 = Kl -JUMP2 

JUMP = JUMP2 
70 K3 = NTLt I ) 

K4 = K3 - JUMP 
KT = OANDWD - K3 
/ L = Kl f KT 

^ A(<3.L) = A(K3*L) +'l.0E*^l5 

. L = K'2 +-KT 

/ A(K3.L) ,= A(K3»L) - l.OE+15 

. I = K4 + KT 

A(K3.L) ^ A(K3,li) +1.0EF15 

A(K3.8AN0WD) = A(K3»BANDW0) - l.OE+15 
80 CONTINUE 
C 

666 FORMAT ( 1 0 X , • At • , I 3 . * . • . I 3. • ) = ’ » E 1 4 . 5 . 1 0 X . * E ( • . I 3 . 
7?7 FORMAT (10E12.4) 

RETURN 

END 


04/25/19 


) = * . E I 4 . 5 ) 


VD 
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SUHROUTINE NUMBER 

THIS SOUBROUTINE RENUMBERS THE NODES AND THE ELEMENTS ON 
SPLIT ING BOUNDARIES SO THAT NODES I. 2 AND 1 ARE ON SPLITTING 
SURFACES* AND REARRANGES THE NODES ON ONE SPLITTING SURFACE 

INTEGER BANDWD 

REAL MSRATE,MSFLUX.LOAD,MACH2 

CaMMQN/ASEMOL/B( A. 4) . VOL{ 3720) .NODE! 4. 372Q) »A( 960,406) ,NTL( 155) , 
PRTIAL(3 ,4 ) , ART ( 155 ), LOCI 155),M0C( 155 ) ,NBC( 155 ) , 

X( 960) ,y( 960) , ZI96 0) .LOAD (960) .BANDWD, E( 1 55 ) , MV ( 1 55 ), 
, NSBf 4, 155), MSI , MS2 . MSI , MSTOT , NTOT , JUMP 1 , JUMP2 , NANZ , 

' NRNZ ,POTEN(960 ) 


no 10 I-1.MS3 
DO 10 J=l,4 
10 N5D( J, I ) = 0 

MS = MSI + 1 

DO 20 J=MS,MS2 
DO 20 I-l ,4 ' 

IF <LBC(J) .EQ. N0DE( I ,MV( J) ) ) NSH(1,J) = NODEII.MVIJ)) 

IF (MBC(J) .EQ. NQDE( I,MV{ J) ) ) NS8(2,J) = NODE ( I , MV ( J ) ) 

IF <N0C(J) .EO. NODE( I ,MV ( J) ) ) NSB(3,J) = NODE ( I , MV ( J ) ) 

20 IF < NODE! I , MV( J) ) .NE. NSB(1,J) .AND. 

. NODE( I ,MV( J) ) .NE. NS8(2,J) .AND. 

. NODE! I ,MV( J) ) .NE. NSB(3,J)) NSB(4,J) = NODE! I , MV ( J ) ) 

K = 0 

MS = MS2 F 1 
DO 40 J=MS.MS3 
DO 30 1=1.4 

IF (LBC(J) .EQ. NODE( I ,MVI J) ) ) NSd(l.J) = NODE ( I . MV ( J ) ) 

IF (MHCIJ) .EQ. NO()E( I ,.MV( J ) ) ) NS8(?,J) = NODE { I , MV ( J ) ) 

IF (NOC(J) .EQ. NOOE( I.MVIJ ) ) ) NSB(3,J) = NODE ( I , M V ( J ) ) 

30 IF ( NODE! I ,MV( J) ) .NE. NSB(l.J) .AND. 

. NlJDE( 1,MV( J) ) .NE. NSB(2,J) .AND. 

. NODF( I ,MV( J) ) .NE. NSB(3,J)) N3B(4,J) = NODE ( I , MV ( J ) ) 

RFAKRANGE THE NODES IN A DECREASING VALUES ARRAY 

00401=1,3 
K = K + 1 

40 NTL ( K ) = NSB( I , J ) 

HSD = (MSS- MS2)*3 
MS01 = MSD - 1 
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NUMBER 


K = 1 

C 

DO 70 L=l *MSDI 
MNP ^ NTL(L) 

00 60 M=L»MSD 

IF (NTLCM) - NTL(D) 60,60.50 
50 K = M 

NTL(L) ^ NTL(M) 

60 CONJINUE 

NTL(K) = MNP 
K - L + I 
70 CONTINUE 
C 

RETURN 

END 

\ 


H 


DAT 


B22 39 


OA/25/19 
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GAUSS 


DATE 


82239 


OA/25/19 


SUBROUTINF GAUSS 

SOLUIION OF LINEAR SYSTEMS OF EQUATIONS RY THE GAUSS ELIMENATION 
METHOD* FOR NON SYMETRIC BANDED SYSTEMS. 

A ; CONTAINS THE SVSTEM MATRIX* STORED ACCORDING TO THE NON 
SYMMETRIC BANDED MATRIX 

POTEN : CONTAINS THE UNKNOWN POTENTIALS AFTER SOLUTION. 

INTEGER 8ANDW0 

PEAL MSRATE*MSFLUX*L0AD*MACH2 

COMMON/ASEMBL/B<4.4)*VOL(3720)»NODE{4.3720) *A(9f»0.406) *NTL( 155) * 
PRTI AL(3*4 ) .ART (155) *LBr ( 155 ) *MBC( 155) .NBC( 155) . 

X( 960) *Y(960) .ZI960) .LOAD (960) .BANDWO.EC 155) . MV ( 155) 
♦ NS0(4.165) .MSl*MS2.MS3.MSTGT.NTfJT, JUMPl. JUMP2.NANZ. 

NRNZ.POTEN(96.0) 


I FORMAT (IHl. 5(/). 10(1H=)**»>> WAf^NING I*./. 

. 15X. ‘SINGULARITY IN ROW*. 15./. 

. ISX. •CALCULATION IS STORED. • ) 


NMI = NTOT - I 
DQ 60 K=l,NMl 
KPl ~ K 1 

IF(ABS( A(K.BANOWD) ) .LT. .00001) GO TO 80 

DIVIDE ROW BY DIAGONAL COEFFICIENT 

NI = KPl + BANDWD - 2 
L = MINOCNI .NTOT) 

DO 30 J^KPl.L 
K2 = BANDWD + J - K 
3 0 A( K »K2 ) = A ( K. K2) /A( K. BANDWD) 

POTFNCK ) = POTLN(K) /A( K. BANDWD ) 

ELIMINATE UNKNOWN FIA(K) FROM ROW I 

DO 50 I=^KP1 .L , 

K 1 = BANDWD + K - I 
DO 40 J=KP1.L 
K2 = BANDWD + J - I 
K3 = BANDWD + J - K 

^0 A(lfK2) = A(I.K2) - A( I .K 1 ) ♦AIK .K3) 

50 f*OTEN(I) - POTFN(I) - AC^ . K1 ) ♦POT EN ( K ) 
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60 CONTINUE 
C 

r COMPUTE LAST UNKNOWN. . 

C ' 1 

IF( ABS( A(NTGTtBANDWO) ) •LT* 0#000000l) GO TO 00 
POTE'N( NTOT) := POTEN( NT OT ) / A ( NTOI f fJANDW D ) 

C 

C APPLY BACKSURSTITUTIOM PROCESS 10 COMPUTE REMAINING UNKNOWNS^ 

. C 

on 70 1=1, NMI 

K = NTOT - I 

KPl = K I 

NI = KPl + BANDWO ~ 2 

L ~ MINO (NI ,NTOT ) 

DO 70 J=KPl ,L 

K2 ~ HANOWD J - K 

70 POTEN(K) = POTEN(K> - A(K . K2 )=<tPOTEN ( J ) 

C 

GO TO 90 
80 WRIT E(6 • 1 > K 
STOP 
C 

90 f^ETURN 
END 
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c 

c 


c 

c 

c 


c 


t 

c 

c 

c 

c 

c 

c 


r 
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SUHROUT INF YSI C 

THIS SUBROUTINE CALCULATES THE PHYSICAL PROPFRTIFS OF THE FLUID 
AT ALL nodes: VELOCITY# VELOCITY COMPONENTS# MACH NUMBER# 

AND NORMALIZED DENSITIES# 

INTEGER BANOWD 

REAL MSRATE# MSFLUX#L0AD#MACH2 ’ , 

COMMON/ASEMBL/BIAtA ) .VOL (37 20) #N0 DE ( 4 , 3720 ) . A ( 96 0 # 406 )\NTL ( I 5S ) # 

• ' PRTIAL(3.4) ,ART( 155 ) #LBC( 155 ) #MnC( 1 55 ) • NBC ( 1 55 ) . 

- X( 960 > • Y(960) #Z( 960 ) »LOAD( 96 0) #BANDWO#E( 155) #MV( I 

• , NSB(4*155) .MSI # MS2 . MS3 # MSTOT • N TO T # JUflP I # JUMP 2 . N AN 

# NRNZ #P0TEN(960) 

COMMON/PHYSI/GAMA.R j ST AGT # ERROR # V ELO ( 3 ) 

C OMMON/INOUT/IND AT A #MTOT .MSRATE.EPS I LN# ITER# I TERMX # NC # ST AGRO # 

# J 1( 60) #NU( 60) # THETA< 60) #ROELEM( 372 0 ) # NN 


EXPO =1#0/(GAMA - 1#0) 

SONIC = GAMA«R*STAGT 
ERROR = 0#0 , 

on 10 1=1 ♦NTOT 
on 10 j=i»4 
10 A( I # J) = 0#0 

DEFINE IHE B MATRIX AND CALCULATE THE VELOCITY COMPONENTS 
VELO(I) (I=l#2#3 FOR X#Y#Z) IN VOLUME ELEMENT II : THE 

potentials ARE NOW IN ARRAY POTEN(J)# IN THE SAME LOOP# 
calculate the element DENSITY AND THE ERROR FOR CONVERGENCE TEST. 

DO 50 II=l#MTOT 

' call V0LELM(2#II) 

V = VOL( I I ) *6 

\ 

DO 20 1=1,3 

20 VELO{ I ) = 0. 0 

DO 30 1=1,4 

on 3 0 J=1 .3 

^0 VELO(J) = COFACTI I , J+1 ,BJ*POTENINODE( I , I I ) )/V f VELO ( J ) 

VSQREO = VELO{ 1)442 + VELO(2)4*2 + VELQ(3)44? 

UENS = (1.0 + VSQRED/I 24SDNIC4EXP0 - VSQPEO ) ) 44 ( -EXPO ) 

ERHOP = ERROR + ( ( ROEL ( I I ) ~ DENS ) /ROELFM ( I I ) ) 4 4 2 

RnELLM(II) = DENS 


C 


\n N 
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on 4 O I = 1 . 4 ^ ' 

AiNODEI I . I I ) f 1 ) = ACNODEI I •! I ) • 1 ) 4^ l.O 
DO 40 J=2*4 

40 AINODEI I f I I ) . J) = A(NODE( I ,I I )♦ J) 4- VHLn(J-l) 

50 CONTINUE 
C ‘ 

IF (ERROR. GT.EPSILN .AND. ITER .UT . ITERMX ) RETURN 
C 

C CALCULATE THE VELOCITY A(If5), THE VC;LaCITY COMPONENTS A(I.J) 

C (J=2.3,4 FOR X*Y.Z)* MACH NUMBER A(1.6), NORMALIZED PRESSURE 

C A(IW)* NORMALIZED DENSITY A{I»8) 

C 

DO 90 I=l,NTOT 
DO <>0 J=2.4 

60 A{ I . JI = A( I . J )/A( I , 1 ) 

C 

VSQRED = A(I*2)*«2 4 A (I . 3 ) ^ 4^2 4 A(I,4)44c? 

A(I,5) = SQRT(VSOREO) 

MACH2 ~ VSQRED/( SONIC -- VS QREO/ ( 24 EXPO ) ) 

A( I , 6) = SQRT(MACH2) 

A ( I ) “ ( 1 #0 t MACH2/ (24EXP0 ) ) 44(-EXP0 ) 

A( I ♦ 7) = A( I t8) 44GAMA 

90 CONTINUE 
C 

C CR I r ICAL VELOCITY 

VrCRI r=SQRT (24STAGT/(EXP04(GAMA41 . 0) ) ) 

PI Z.3. 1 41 6 
C 

r CALCULATE THE RADIAL AND THE TANGENTIAL VELOCITIES 

C ' TANGENTIAL VELOCITY IS STORED IN A(I,2) AND THE F^AOIAL 
C VELOCITY STAoRD IN A(It3) 

C 

DO 1 00 K=l • NC 
JK-J 1 { K) 

JKI=J1 (K)4NU(K)-1 
DO 100 I=JK,JKl 

VELD XY=SQRT( A( I • 2) 4424A( I , 3) +4 2) 

PHI=;ATAN(A(I.3)/A( lp2\) 
f>SI = THETA( K)--PHl 

IF ( A( I,2)4A( If 3) .LT.0.0) PSI=P1-PSI 
A { 1 .3 )=~ ( A( I f 3 )/AHS ( A( I * 3) ) ) 4V ELOX Y 4COS ( PS T ) 

100 A(rf2)=~(A(T*2I /AOS( A( 1*2))) 4VFL0XY4SI N( PS I ) 


C 


Rh TURN 
f^ND 
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FUNCTION C0FACT( If J»B> 

FUNCTION COFACT CALCULATES THE COFACTORS OF THE MATRIX B(I,J) 
WHERE dCI.D-l ANQ J IS ANY REAL NUMBER 

DIMENSION B(4t4)» C(4,4) 

DO I M = 1 .4 
00 I K=U 4 

1 CCK,M) = B(KfM) 

IF ( I .EQ. 4 ) GO TO 3 

REDUCE 4t4 MATRIX TO 3*4 MATRIX 

DO ? M=1 ,4 
DO 2 K=I *3 

2 C(K^M) = CC (KFl ) .M) 

3 IF ( J .EQ* 4) GO TO 5 

REDUCE 3+4 MATRIX TO* 3^3 MATRIX 

DO 4 M=J*3 
DO 4 K=l *3 

4 C(Kf M) = C(K, (MFl) ) 
evaluate the COFACTOR 

5 DETERM = C ( I , 3 ) * ( C ( 3 • 2 ) - C<2,2)) F 

• C (2%3) F(C( I ,2) - C(3.2)J F 

• C(3»3)*(C(2»2) - C(lt2)I 

L = (-1 ( IFJ) 

rOFACT = DETERM F L 

RETUF?N 

END 


L 
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SUBROUTINE VOLELMIN.K) 

THIS SUBROUTINE ElUILDS UP THE MATRIX OII.J) WHf. RE Q(I.l) = l.O 
N = 0. CALCUATES THE VOLUME VOL(K) OF AN ELEMENT K. 

N = 1, UUIL UP MATRIX B FOR AN 'ELEMENT WITH SPLITING BOUNDARY 

N = 2. BUILD UP MATRIX 0 FOR A VOLUME ELEMENT 

INTEGER BANDWD 

REAL MSRATEt MSFLUX .L0AD.MACH2 

C0MM0N/ASEMBL/B(4i4) . VOL { 37 20 ) , NODE ( 4 t 372 0 ) t A { 96 0 . 406 ) t N TL ( 155), 
PRTIALI 3,4) ,ARTC 155) ,LQC( 1 55) . MBC( 155) ,N3C( 155 ) , 

X (960 ) , Y <960 ) , Z( 960 ) .LOAD ( 960 ) . BANDWD, E( 155) ,MV( 155) 
. NSB(4,155) ,MS1 ,MS2 , MS3 ,MSTOT ,NTOT , JUMPl , JUMP2.NANZ, 

NRN2,P0TEN(960) 


IF ( N ,NE. 1 ) GO TO 2 
DO 1 1=1,4' 

B ( I , 1 ) = 1.0 
B( I , 2) = X( NS«{ I ,K ) ) 

B( I ,3) = Y(NSB( I,K) ) 

1 B( I , 4) , = Z( NSB{ I ,K) ) 

RETURN 

2 DO 3 1=1,4 

B( 1 , 1 ) = 1,0 

8(1,2) = X(NOD£( I.K) ) 

B( I , 3) = Y( NODE( I ,K) ) 

3 8(1,4) = Z ( NODE ( I.K) ) 

IF (N ,EQ. 2) RETURN 

C 

VOL(K) = 0.0 

DO 4 1=1,4 

VOLL = COFACT( I ,2,0)«R( I ,2) 

4 VOL(K) = VOL(K) + VOLL 
VQL(K) = A8S(VOL(K ) )/6 

C 

RETURN 

END 
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A a 


FUNCTIGlsr AR( I ) 

THIS FUNCTION CALCULATES THE AREA OP A TRIANGLE IN CARTESIAN 
COORDINATES. l = NUMDER OF SURFACE ELEMENT UrvlUFR CONSIDERATION. 


INTEGER BANOWO 

REAL MSRATEf MSFLUX,LOAD.MACH2 

CnMM0N/ASEMBL/8(4»4) tV0L(3720) . NODE ( 4 . 3720 ) .A (960.406 ) .NTL( 155) . 

PRTIAL(3,4)» ART(155).LBC( 155) . MRC ( 155) .r4BCC155) . 
X(960)»Y(960),2C960) *LOAD ( 960 ) . FMNDWD. E( 155) .MV( 155) 
. N5BC 4, 1 55) .MSI .MS2. MSI . MSTOT .NTOT. JUMP I . JUMP2.NANZ. 

NRN2.P0TEN(960) 


N) 

CO 


Al 


A2 


A3 

AR 


Y(LBC( I ) >»(Z(MI3C{ n ) 
Y(MHC( I ) )♦( Z( NBC( I ) ) 
Y(NRC( I n*(Z(l OC( 1 ) ) 
X{LBC( I ) )»(7(MBC( I ) ) 
X(MBC( I ) l*(Z(NBC( I ) > 
X(NQC( I ) )*(Z(LBC( 1 ) ) 
X(LHC( I ) )#(Y {MBC( I ) ) 
X(MBC( I ) )<<( Y( NBC( I ) ) 
X( NBC ( I ) )♦( Y(LOC( I ) ) 
SQnT(Al**2 + A2**2 +■ 


- Z(NBC( I ) ) ) + 

- Z(LBC( I ) > ) + 

- Z( MBC< I ) ) ) 

- 2{NBC(n>) + 

- Z(LBC( I ) ) ) + 

- Z(MOC(I)>) 

- Y(NBC{ I ) ) ) + 

- Y(LBC{ I ) ) ) 4- 

- Y(M0C{ I ) n 
A3*»2)/2 


C 


RfTURN 

FND 



i 
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K> 

KO 


C 

c 

c 


c 


c 

c 

c 


SUBPnUTlNE INPUT 

THIS SUDPOUTINE READS AND WRITES THE INPUT DATA. 

INTEGER BANDWO 

REAL MSRATE,MSFLUX*LOADf MACH2 

COMMON/ASEMBL/B( A #4 ) , VOL (372 0 ) • NODE (4» 37 20 ) . A ( 96 0 , 4 06 ) • N TL < 1 55 ) ♦ 

. PRTI AL(3t 4) • ART( I 55) *LRC( I 55) • MUC (155) *N0C( 155) » 

. ' X (960 ) ♦ Y(960 ) • Z(960 ) • LOAD ( 96 0 ) * OANDWO.E( 155) • MV ( 155) 

. f NSB(4.155) »MS1 * MS2 . MS3 ♦ MST OT . NTOT . JUMP I . JUMP2.NANZ. 

. NRNZ.POTEN(960) 

COMMON/INOUT/INOATAtMTOTtMSRATE.FPS ILN, ITER» I TEPMX ♦ NC * ST AGRO • 

. J1 (60) .NU(60) »THETA(60) . ROELEM ( J72 0 ) . NN 

COMMON/PHYSI /GAMA, Rt STAGT TERROR .VELOI 3) 


1 FORMAT 

2 FORMAT 

3 FORMAT 
a FORMAT 

9* FORMAT 


( 6F 10. 0,215) 

(1415) 

( 1P5E14.5) 

( IHl ,////, SOX, • I N P U T 

SOX, • 

( 33 X,* TOTAL NUMBER 0P“ 
33X, MOTAL NUMBER OF 
33X, • 

33X , • 

3-^X, •SURFACE 
33X,* 

33X , • 

33X, • 


DATA*,/, 

! • ,////) 

elements'^Tmtot ) : • , 

NODES (NTOT) :» , 


ELEMENTS AT BOUNDARIES 


( 

( 

< 

( 

( 

( 


33X, • BANDWIDTH OF 5T1FNESS MATRIX 


, 16,/ , 

, 16 •/// 
MSI :• 
MS2 : • 
MS 3 :• 

mstot: • 
nodin: • 

JUMP2 : • 

* - 


33X, ‘MAXIMUM NUMBER OF ITERATIONS 
33X,«TOTAL NUMBER OF SFCTIONS 
33X, • SPECIFIC HEAT RATIO 
33X, ‘CONVERGENCE CPITFRION 
33X, • STAGNAT ION TEMPERATURE 
33 X.* STAGNATION PRESSURE 
33X, ‘IDEAL GAS CNSTANT 
»// ,33X , • MASS PATE THROUGH SYSTEM 


( BANDWO) 

{ I TEPMX) 

( NC ) : • ♦ 

(GAMA) : ‘ , 

(EPS ILN) 

( STAGT) : • ,F1 1 
(STAGP) : ‘ ,Fl 1 .5 
( R) : • ,F1 1 .5 

(MSRAVr ) : • ,FI 


,15,/, 

, 15,/, 

•15,/, 

,15,/, 

• 15,/, 

.15,///, 

15,//, 

IS,//. 

15.//, 

FI 1.5,//, 
F12.6,//, 

.5, • R • , //, 
PSF‘,//, 
FT2/S2.R‘ 
I . 5, • SLG/5‘ ) 


RFAD( 5 , I ) 
RFA0(5,2) 
RrAO(5,2 ) 
PFAD(5,2) 
M = MSI + 
READ (5 ,2 ) 
RFADC 5. 2 ) 
READ (5,2) 


M SR ATE ,R .GAMA, STAG Tf SI AGP.E'^SILN. INDATA, ITERMX 
MTOT, NFOT .MSI ,MS2, MS3. MSTOT .NAHZ .NR N/ . NC .BANDWO 
(LHC( I ),MBC( I ),NBC( I ) , I=l .MSrOT) 

( NUDE( 1 . J> .NODE ( 2, J) .NODE ( 3 . J ) .NODE . J ) , J- 1 ,MTOT ) 
I 

(MV( I ) , l=M.MS3 ) 

( J 1 ( I ) . I =n . NC ) 

(NU( I ) , 1=1 ,NC> 



nnnnn n onon 


AN IV G1 RHLHASE i> , 0 


INPUT 


DATE 


8?219 


04/25/19 


RfAD( 5 ,3 ) 
READ (5.3 ) 
PEAD(5.3) 
REAIX 5, 3) 


( THETA( I ) . 1=1 ,NC) 
(X( I) , I=l,NTOT) 

CY ( n , 1=1 .NTOT ) 
(Z( I) .1=1 .NTOT) 


INITIALIZE NORMALIZED ELEMENT DENSITIES TO HE UNITY AND CALCULATE 
THE STAGNATION DENSITY. 

DO 50 1=1. MTOT 

50 ROfLEMI I ) = 1.0 

STAGRO = STAGP/I R*STAGT ) 

JUMPl = NANZit'NRNZ 

JUMP 2 = LBC(MS2+n - LBC(MSl + n 

IF (BANDWD .LT’. JUMP2) BANDWD = JUMP2 

IF < INDATA .EQ. 0) RETURN 


WRITE INPUT DATA IF FLAG INDATA=l 
WRITE(6.0) 

WRIT E( 6.9 ) MTOT. NTOT .MS I , MS2 . MS3 . MSTOT , NU ( 1 ) . JUMP2 . 0 ANDWD . ITFRMX. 
NC .GAMA.EPSILN.STAGT.S1 AGP.R.MSRATE 


129 RETURN 
F ND 


C 



U o - u u 


AN IV G1 RELFASE 'i.O 


OUTPUT 


DATE 


82239 


04/25/19 


sunnouTiNE output 

c 

C THIS SUBROUTINE PRINTS OUT THE CALCULATED PHYSICAL PROPERTIES 

C AT THE END OF EACH ITERATION* 

r 

INTEGER BANDWD 

REAL MSRATE , MSFLUX .LOAO.MACH2 

COMMON/ASEMQL/B( 4.4) , VOLC 3720) .NODE { 4.3720) .A( 960,406) .NTLU55) . 

. PRT I ALI3.4 ) . ART ( 155 ).LBC( 1 55 ) . MDC ( I 55 ) .NBC{ 155) , 

. XI 960) , Y( 960) , Z(960) .LOAD! 960 ) .BANOWD, E( 155 ) ,MV( 155) 

. , NSB( 4,155) .MSI ,MS2, MS3. MSTOT.NTOT, JUMPl . JUMP2 »NANZ, 

. NRN7.POJEN(960) 

, COMMON/INOUT/INDATA.MTOT.MSRATE.EPSILN, I TER , I TERMX , NC , ST AGRO , 

. J II 60) ,NU( 60) , THETA I 60) .ROELEMI 3720) , NN 

COMMON/PHYSI /GAM A , R , ST AGT , ERROR , V FLO ( 3 ) 


1 FORMAT (1415) 

2 FORMAT I IHl ,//,7X, NODES POTENTIALS VELOCITY-T VEL*.» 

•OCITY-R VELOCITY-Z VELOCITY MACH NUM. NORM. PRES.*. 

. * NORM. DENS. • ,/) 

3 FORMAT I 19, 17, 1P8EI4.S) 

4 FORMAT (IP5E14.5) 

5 FORMAT (2X,*END OF SECTION*,I3) 

6 FORMAT (8(/),43X.*E ND OF PROGRA M • , /. 48X , 2 71 I H_) . 

. ////.33X, ’CALCULATIONS ENDED :*, 

. //, 33X, • NUMBERS OF ITERATIONS PERFORMED : ITER =*,I3, 

. // ,33X , ’COVFHGENCE REACHED : ERROR =*,lPEa.l) 


’WRITE! 6, 6) ITER, ERROR 
J = 0 

L = 53 

DO 9 K=l,NC 
JK = J 1 (K ) 

NK = J1 IK) + NU(K) -I 
C 

DO 8 I=JK,NK 
J = J + I 

L = L + I 

IF (L - 5J ) U.8. 7 

7 L = I 

V'» P n E ( 6 , 2 ) 

3 ^a.R1TE(6.3) j. I .POTENI I) , I A( I ,M) . M=?,0) 
ARirE(6,5) K 
• 9 CONI INUE 



XN IV Gl Rfc'LEAfilf ?«0 


OUTPUT 


DATE = 82239 


C 

C 

C 

C 

C 


C 


WRITE(7 .1) 
WRITE! 7» I ) 
WRITE{7,4) 
WRITE! 7 .4 ) 
WRITE! 7,4 ) 
WRITE!7,4) 
WR I T E ! 7 , 4 ) 
WR IT E! 7 ,4 ) 
WRITF!7,4) 
WRITE! 7,4 ) 
WRITE! 7,4) 


(Jill), 1= 

<NU!I), 1 = 

(X(I), I=l 
(Y(I), 1=1 

( Z( I ) . 1=1 

( POTENC I ) . 


( A( I ,6) 
( A{ 1,7) 
!A(I,2) 
( A! I ,3) 
( A! 1,4) 


1 ,NC) 

1 ,NC ) 

,NTOT ) 
,NTOT) 
,NTOT) 

1=1 .NTOT ) 
= 1 .NTOT) 
=l,NTOT) 

=1 ,NTOT ) 


I 
1 
I 

1 = 1 ,NTOT) 
I 




RETURN 

END 


/ 


OJ 

K> 


04/25/19 



NUMERICAL EXAMPLE 


A numerical example is presented to illustrate the use of 
the program and show the results of the computations. The 
flow IS investigated in a radial inflow turbine for which the 
experimental measurements were reported in reference [8] . The 
mass flow rate is one Ib/sec (0.454 kg/sec) and the scroll 
inlet section diameter is 4.88 inches (1.92 cm). The niomber 
of nodal points used in the discretization pattern was 956 
placed on 24 scroll cross sections as shown in Figs. 3 and 4. 

The number of tetrahedral elements resulting from this dis- 
cretization pattern was 2,716. The printed output from this 
example is given in Table 2. The first part of the output 
consists of the performance parameters and other input data. 

The second part of the output represents the results of the 
computations and consists of the velocity potential, the velo- 
city components, the Mach number, the normalized density and 
pressure. The pressure and density are normalized with respect 
to their inlet stagnation values. The output data from the 
computed results in the various scroll cross sections. The 
contours of the computed velocity potentials are shown in 
Fig. 5. The jump in the value of the velocity potential across 
the splitting surfaces was found to be equal to 150 as can be 
seen in Fig. 5. The same figure illustrates the same slope 
for the potential contours on the opposite sides of the splitting 
surface and hence the continuation of velocity pressure and 
density fields. The contours of the computed velocity component 
in the circumferential direction are shown in Fig. 6. The 
measured [5] circumferential velocity components are shown in 
Fig. 7. In both Figs. 6 and 7, the values of the velocity 
componet are normalized with respect to the inlet velocity. 

The disagreement between the computed and measured results 
is mainly in the boundary layer region over the outer scroll 
surfaces. One can also observe some secondary flow development 
in the measured results. However, both boundary layer and 
secondary flow effects decrease with increased velocities. 
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The agreement between the computed and measured results near 
the scroll neck, in spite of the very low velocities in the 
case presented, with an inlet velocity of 41 m/sec, leads to 
the conclusion that the velocity and pressure fields in this 
region are not very sensitive to these effects. The flow is 
greatly accelerated in this region and it also acquires a 
radial velocity component as evidenced by the shape and density 
of the potential lines in the vaneless nozzle portion in 
Fig. 5. 
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TABLE 1 


INPUT DATA FORMAT 


Type of Input Data 

Input Variables 


Format 

Performance, Printout and 
Convergence Parameters 

MS RATE, R, GAMA, STAGT, 
INDATA, ITERMX 

STAGP, EPSILON, 

6F10.0, 215 

Discretization Parameters 

MTOT, NTOT, MSI, MS2, MS3, MSTOT, 
NANZ, NRNZ, NC, BANDWD 

1415 

Topology of Permanent 

Surface Elements 

(Inlet, Splitting, and Exit) 

(LBC(I), MBC(I), NBC(I) 

, 1=1, MSTOT) 

1415 

Topology of Volume Elements 

(N0DE(1,J), N0DE(2,J), 
J=l, MTOT) 

N0DE(3,J), N0DE(4,J), 

1415 

Volume Elements Affected by 
Splitting Surfaces 

(MV(I), I = MSl+1, MS3) 


1415 

Cross-Sectional Data 

(J1(I) , I = 1,NC) 
(NU(I) , I = 1,NC) 
(THETA (I), I = 1,NC) 


1415 

1415 

1P5E14 . 5 

Nodal Coordinates 

(X(I) , I = 1, NTOT) 

(Y(I) , I = 1, NTOT) 

(Z(I) , I = 1, NTOT) 


IP5E14.5 

It 

If 



INPUT DATA 


TOTAL NUMBER OF ELEMENTS (^MTOT) : 3716 
TOTAL NUMBER OF NODES (NTPT) : 956 



{ MSI 

48 


( MS 2 

64 

SURFACE ELEMENTS AT BOUNDARIE 

S < MS3 

30 


( MSTOT 

152 


{ NOD IN 

35 


( JUMP2 

19 

BANDWIDTH OF STIFNESS MATRIX 

(BAMDWO ) 

203 

MAXIMUM NUMBER OF ITERATIONS 

( ITEPMX) 

5 

TOTAL number OF SECTIONS 

(NO 

24 

SPECIFIC HEAT RATIO 

( GAMA ) 

1 .40000 

CONVERGENCE CRITERION 

( EPSILN) 

0^.000100 

STAGNATION TEMPERATURE 

(STAGT ) 

530.00000 R 

STAGNATION PRESSURE 

(STAGP) 

4319.67969 °SF 

IDEAL GAS CNSTANT 

(R) 

1716.26001 FT2/S2.R 

MASS RATE THROUGH SYSTEM 

(MSRATE ) 

0.03105SLG/S 


END OF PROGRAM 
CALCULATIONS ENDED : 

NUMBERS OF ITERATIONS PERFORMED : ITER = 2 

COVERGENCE PEACHED : ERROR = 5.3E-06 


TABLE 2: OUTPUT DATA, FIRST SET. 
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ru)Of r> 


u> 

03 


/t 

1 1 

I 

J *3 

^ 4 

n r» 

6 

7 7 

It 


9 

9 

1 0 

10 

1 1 

1 1 

12 

12 

1 3 

13 

14 

14 

1 5 

1 5 

15 

16 

1 7 

IV 

18 

18 

19 

19 

20 

2 0 

?l 

21 

2? 

?2 

?J 

23 

24 

24 

?5 

25 

?5 

2f| 

27 

27 

28 

28 

?9 

?9 

JO 

30 

31 

31 

32 

32 

33 

33 

34 

34 

35 

35 

Of Sf 

CT ION 

Jt> 

36 

^ 7 . 

37 

3 8 

38 

39 

.39 

40 

4 0 

4 1 

4 1 

4 ^ 

4 2 

4 J 

4 i 

44 

4 4 

4 5 

45 

46 

46 

47 

4 7 

48 

4 8 

4 9 

49 

50 

50 

51 

5 1 

5? 

52 

5 3 

5 3 


IM)Trr4T I ALS 

“ I . ooooor40t) 

- I • OOOOOf 400 
~l . 000001*400 
“I • OOOOOf 400 

- I . OOOOOf 4 00 
-I .OOOOOf 400 

- 1 • onooor ton 
- 1 * ooooor 400 

-1 .OOOOOE-400 
- l . OOOOOf 4 00 

- I . OOOOOf 400 
- 1 • 000001 400 

- 1 • ooooof'f no 
-I .ooooor 4 00 
- 1 « noooorfoo 

- I .OOOOOt 400 

- l • OOOOOf- 400 
-1. OOOOOf 400 
-> I • OOOOOt 4 00 
-I.OOOOOF400 

- 1 • ooooor 400 

-l.OOOOOLfOO 
-I . ooooor.400 
- 1 . 000001:400 
~l .OOOOOF400 
-1 .0000UF4-O0 
-I.OOOOOF400 
-I .OOOOOF400 

-i.ooooor+oo 

-I.OOOOOFfOO 

-i.ooooor40o 

-1. 00OO0F400 

-i.ooooor+oo 

-U OOOOOf 400 

- 1 . ooooor 400 

1 

7. U)7b3E400 
7, l*^864F+00 
7. I 7?a 1F400 
7. I6706E+00 
7. 16065F400 
7. IS?OTL400 
7. I 3969E400 
7. I 5037r+00 
7. 14420L400 
7* I2620r 400 
7. 1 1 2894 +00 
7. I 0 'j 96L 400 
7.1 0050F 400 
7. 10890F400 
7. I 231 7F 4 00 
7. I l230CfOO 
7. 09f476E400 
7. 0044 9f 400 


vFi nr iTY-r 

4 • l40Sir 401 
4 . 141 I 7f.4 0l 
4 . I i04)4f 401 
4 . Ml *-.2t fOl 
4 • l?dl 3F401 
4. l?S0Hf:401 
4 . 12419F + 01 
4. 13272F + 0I 
4. ir73lr +01 
4. I 1929F+01 
4. ll?94r+01 
4 . ionopL + 01 
4.10669rf0l 
A.llO^^OFtni 
4.11 103E + 0I 
4. 1103SE40I 
4. 1034 7F+01 
4.096^»3F + 01 
4.0«0?‘SEf 01 
4.005f3L+01 
4 .004?2E+01 
4.08 75SE401 
4.0H559E+01 
A .07929r+01 
4.0730JE+01 
4.0ftfif>3C40l 
4.00521E+01 
4.06306C+01 
4 .07004H401 
4 • 06371F: + 01 
4.05*>44f: + 01 
A *05079F+0 I 
4.0500JE40I 
4.05164E+0I 
4 .0S270F+01 

4.09095E+01 
4.07953E+0I 
4 .0536*»E + 01 
4.04 I 20E40I 
4 ,03905Ef01 
4 .04340E+01 
4.0fi061E + 01 
4. I3231b+01 
4.1 2340E+01 
4. I0419E + 01 
4.0924 lE+Ol 
4 .00363Ef 01 
4.O7710E+0I 
4 .O075OE+O 1 
4.1 7365E+01 
4.l4fl92C+Ql 
4.13442F+0I 
4 . I 1905E + 01 


VELftri TY-fl 

-2.77710E-0I 
.3?lfi7r-0l 
-2. b«ft2HL-0l 
-2.015ML-01 
-3.0M6IF-01 
-3.90991F-0I 
-*5.J03r>0E-0l 
-I . 7509Me-0l 
-I • 10277E-01 
-9.12454E-02 
-l.03‘539F-0l 
-1-34274F-01 
-1.091 lir-01 

- I .6659IE-01 
-I .94307r-01 
-1. 38753r-0l 
-1 .22124E-01 
-1 .19I43E-01 
-1 .40290E-01 
-1.7ri4b7F-0| 
-I .407S4E-01 
-3. 305??e-01 
-2.6731RE-01 
-2.43147E-01 
-2.2236«e-0l 
-1 .99H09E-01 
-1.98640E-01 
-1.44093E~01 

1 .77200F-01 
1 • 29026E-04 
-1.03290E-04 
1.27106F-O5 
-1.03I52F-04 
1.272136-05 
1.266776-04 

-3.96642F-02 
1 « 04591E-0I 
3. 1 79016-02 
-I *632826-02 
-1. 100386-01 
-3.559136-01 
-5. 70675C-01 
-1.076266-01 
-8.288346-0? 
-3.486146-02 
-0. 1 40026-02 
-I .700986-01 
-3.269946-01 
-5.044046-0 1 
-2.856546-02 

- 1. 120426-01 
-7.300946-02 
-7. 094456-02 



VELOCITY 

-/ 

-2 

• ?8283c“ 

01 

-2 

. 112856- 

0 1 

-1 

. 1 8984*^- 

01 

-2 

• 905716- 

0 f 

7 

.74729F- 

02 

2 

• 2 154 06- 

01 

4 

.27450E- 

0 1 

-9 

• 164?9F- 

0? 

-1 

.013236- 

01 

-8 

.28440E- 

02 

-5 

• 57j781 F- 

0? 

-1 

. 1 54 71E- 

0? 

1 

• fl7591F- 

02 

2 

. 284056- 

01 

-? 

.576456- 

0? 

-6 

.066 7 0*^- 

0^ 

-6 

.070556- 

0? 

-7 

.684493- 

02 

-6 

.07,8576- 

02 

-5 

.290096- 

0'> 

? 

•5749?F- 

0 5 

-2 

• 64788E- 

03 

-1 

.3??09r- 

02 

-3 

.015656- 

02 

-5 

. 7.3662E- 

0? 

-7 

•72312E- 

0? 

-7 

.894076- 

07 

3 

• 52859E- 

05 

-3 

.53218F- 

01 

-7 

.391 17F- 

02 

-4 

. 57,4 09F- 

02 

-1 

•29428E- 

02 

1 

.379166- 

02 

1 

•68006E- 

02 

7 

.34329F- 

05 

-4 

.94881F- 

01 

-5 

. 02956F- 

01 

-? 

.801 136- 

01 

-2 

*004876- 

0 1 

-1 

.380976- 

03 

1 

*417636- 

01 

2 

.314416- 

01 

-3 

. 804046- 

01 

-4 

« I9080E- 

0 1 

-3 

.296226- 

01 

-? 

• 42535K- 

01 

-1 

.000026- 

a 1 

-4 

• 34 8 07>E- 

02 

1 

.912I7F- 

01 

-4 

« 027866- 

01 

-2 

.673696- 

01 

-3 

.340006- 

0 1 

-3 

.185916- 

01 


vn UCITY 

4.1 45666+01 
4.141286+0! 

4 . I 15 746 + 01 
4.131626+01 
4 .123256 + 01 
4.125096+01 
4. 124756+01 
4 .132776 + 01 
4.127336+01 
4*1 19116 + 01 
4 • 1 M966 + 0I 
4. 108846+01 
4 .105736 + 01 
4.1 10536+01 
4.1 13086+01 
4.1 10376 + 01 
4.1 03506+01 
4.096556+01 
4.090236+01 
4.085376+01 
4.08425E+01 
4 .087696+01 
4.085656+01 
4 .078376 + 01 
4.073106+01 
4.063596+01 
A .065276+01 
4.063096+01 
4.07103E+01 
4.063726+01 
4.055446+01 
4 .050796 + 01 
4.050036+01 ' 
4.05164C+01 
4 .052706 + 01 

4.099266+01 
4 .079356 + 01 
4.053756 + 0 I 
4.04125F+01 
4 .039076+01 
4.043586+01 
4.069006+01 
4. 133006+01 
4.123706+01 
4.104326+01 
4.092496+01 
4.083686+01 
4.077246+01 
4.087946+01 
4 .173846 + 01 
4.149026+01 
4.134366+01 
4 .1 191B6+0I 


MACfl NUM. 

3.675066-02 
3.07029E-02 
3.665386-02 
3.061726-02 
3.650736-02 
3.656826-02 
3.655636-02 
3.66274C-0? 
3.657926-02 
3.05081E-02 
3*645186-02 
J.64153E-02 
3.639666-02 
3.64303E-02 
3.64599F-02 
3.642896-02 
3.6J679E-02 
3-.630636-02 
3.625066-02 
3.620726-02 
3.619736-02 
3.622776-02 
3. 620996-02 
3.614516-02 
3.609846-02 
3.605056-02 
3.602906-02 
3. 600976-02 
3.600016-02 
3.601536-02 
3.594196-02 
3.590076-02 
3.50939E-02 
3.590826-02 
3.591766-02 

3.6330 36-02 
3. 615836-02 
3.592696-02 
3.581616-02 
3. 579686-02 
3.533676-02 
3.606206-0? 
3. 662956-02 
3.654706-02 
3.637526-02 
3.627036-02 
3.619226-02 
3.613516-02 
3.623006-02 
3.699156-0? 
3.677156-02 
3. 664336-02 
3.650696'-02 


NURM. PRES. 

9.99056E-01 
9. 99059C-OI 
9.99063E-01 
9.99063E.-01 
9. 990G6F-01 
9.99066E-01 
9.99066F-01 
9.99063E-01 
9.99066C-01 
9.99069E-01 
9.99073r-01 
9.99073E-01 
9.99076E-01 
9.99073E-0I 
9.99073E-01 
9.99073E-0I 
9.99076E-01 
9.990796-01 
9.990B3C-01 
9.990866-01 
9.990O6E-0I 
9.99083E-01 
9.99086E-01 
9.99009E-01 
9.99089E-01 
9.99093E-01 
9. 990936-01 
9.99096E-01 
9.99089E-01 
9.99093F-01 

9.990996-01 
9.99099E-01 

9.990996-01 
9.99099E-01 

9. 990996-01 

9.99079E-01 
9.990866-01 
9.900996-01 
9.99103C-01 
9.991066-01 
9.991036-01 
9.99093C-01 
9.99063E-01 
9.990666-01 
9. 990766-01 
9.990836-01 
9. 990666-01 
9.99089E-01 
9.990856-01 
9. 990466-01 
9.99056E-0 I 
9.990636-01 
9.990696-01 


NORM. DENS. 

9. 993266-01 
9. 993206-01 
9.99330E-01 
9.99330F-0I 
9.993336-01 
9.99333E-01 
9.99J33E-01 
9.99330E-01 
9.99333E-01 
9.99335E-01 
9,Q933rE-01 
9.99337E-0I 
9.99340E-0I 

9.993376- 0 I 

9.993376- 01 

9.993376-01 
9.99340E-0I 
9.99342H-01 
9.993456-01 

9.993476-01 

9.993476-01 
9.993456-01 

9.993476-01 

9.993496-01 

9.993496-01 
9.993S2E-01 

9.993526-01 
9.99354E-0I 
9.99349E-01 

9.993526-01 
9 .993576-01 

9.993576-01 

9.993576-01 
9 .993576-01 

9.993576-01 

9.993426-01 
9.99347E-01 

9. 993576-01 
9.99359E-01 
9.99361E-01 
9.993596-01 

9.993526-01 
9.993306-0 I 
9.9933 ^E-0l 
9.993406-01 
9 .R9345H-0 I 

9. 993476-01 

9.993496-01 
9.99345E-0I 
9.993186-01 
9.993266-01 
9.993^06-01 
9.99335F-01 


OUTPUT DATA, SECOND SET 



ff S niULMI lALS yLLUCIl/-r Vr-.LCJC I TY~K 

S/| 5 '4 7.0 71 00 4.l0 3?tCf01 - I . 4 7^> I I 

SS SS 7.or>?70r +00 4.0MO4 3r+0l 6 7'>46t:-0 l 

Sf' *M> 7. ’,r +00 4 .()f‘l?ic + 01 -3 I 7e-0 I 

*■^7 S7 7.07337r+00 4.'’03fllF+01 -3.e37ll3F-Ol 

‘>0 7. 07S3f>f:*00 ) . 104'iOF+Ol -4.003b5r-01 

; *>0 7. Ob^jnot +00 4.14O41C + 0I “4.b6743E“0l 

00 f>0 7.040571+00 4.IJ140C+01 -4.24704F-01 

M Ol 7.03'307F+00 4. 1 1347E + 01 -3.65^9lF-0| 

7.^ 6? 7.0P3 70F + 00 4. 0037 7i: + 01 -3.2B710C-0 1 

r>3 63 7.00S410F.+00 4.07706E + 01 -3.47I45G-0I 

74 64 7.024427+00 4.17335F+01 -6.340lbE-0l 

06 6b 7.007127+00 4.13B02E+01 -6.08900E-01 

<6 66 6.9090Ot +00 4.12168E+01 -6.226077-01 

67 67 6.O770bF+0O 4 . I 0'5 7s>E + 0 1 -4.74616E-01 

63 6U 6.974007+00 4.09390E+01 -4. 367967-01 

60 69 6.978197+00 4.082687+01 -3.86306E-01 

70 70 6.960417+00 4.07109E+01 -3.84171E-01 

END Of 61 CT ION 2 

71 7l l.r.^l‘>2r^ 01 3.912447 + 01 3.18266C + U0 

72 72 1.610317+01 3.88I90F+01 3.389606+00 

/3 7i 1.499777+01 3.83260E+01 ' 3.24066C+00 

74 74 1.494-SOL + Ol 3.61 + 37E + 01 3. 156426+00 

76 76 1.4 9 31 Of +01 3.826907+01 2.939117 + 00 

76 76 1.496377+01 3 . 8+>2 1 f.E + 0 1 2.260467+00 

77 77 l.60Ot»9F + 0l 3.90766F + 01 1 .68157E + 00 

78 78 l.52496r+0l 4.031767+01 2.47146E+00 

7‘> 79 l.620i>47+0l 4.0293 »7+0l 1.982baC+00 

oO 80 1.610P4E+01 3.991J2F+01 2.10123E+00 

'1 81 1.605767+01 3.97281E+01 , 1. 980987+00 

62 02 1.602227+01 3.96620C+01 1.70I62C+00 

P3 03 1.602677+01 3.O7274E+01 1. 306997+00 

04 04 1.601037+01 3.97716E+01 2.97724E-01 

86 85 1.63616F+01 4.19397E+01 2.927927+00 

86 06 1.529777+01 4.I7J92E+01 2*01l43E+00 

, 87 87 1.624067+01 4. 160647+01 2.014867+00 

^ 60 08 1.618407+01 4. I 4 0627+ 01 1.9U>llE+00 

09 09 1.51 POOL +0 1 4.1I694E+01 1.69564E + 00 

90 90 • 1.60 8327 + 0 1 4. 093007+0 1 1.3J544F.+ 00 

91 91 1.60044++01 4 . 06759'>E+01 4.902477-01 

92 92 1.64681F+01 4.37607E+01 2. 397327+00 

9 3 9 3 1.63 80 01+0 1 4.33*>3lE + 0l 1.36916E + 00 

94 94 I.63393F+DI 4.I3624E+01 1.09032C+00 

99 95 I .5271 ?i:+0 1 4.JI0l0t+0l 1.015487+00 

9f 96 l.6?040f+01 4.266607+01 1.05694E+00 

‘^7 97 l.5in0f+01 4.2O046L+OI 1.09141E+00 

98 98 1. 602671 *01 4.1367IE+0I 6.64901F-0I 

99 99 1.646441+01 4.4267»rC+0| 4.040737-01 

K»0 lOO 1.639641+01 4.1745OE+01 1.662617-03 

lOl 101 1.6J045F+01 4. 57»860E + 01 -4.36764E-01 

lt'2 102 1.521927+01 4.33793E+0! -5.61985E-01 

lOI 103 1.6I660E+0I 4.293707+01 -5.1B822E-01 

104 104 1.611077+01 4.23267E+01 -2.90707E-0! 

106 106 1.606507+01 4.16o46r+0l -1.71164F-01 

END nr • srcT I on i 

18'^* 106 2.66443F+01 3*691 14F + 01 3. 019397 + 00 


VELOCITY-/ VILOCITY M4CM MUM. NOPM. PRFS. NORM. DENS. 

-3.0 19l4F-tJl 4.103367 + 01 3.63666L-0? 9.99076E-01 9.993407-01 

662006-01 4.086607+01 3.62161E-02 9.99O03E-O1 9.993457-01 

-2.60 W>2r-'12 4.081 37C+ 01 3.617iaE-0? 9.990067-01 9.993477-01 

-4.426987-01 4.20422E+01 3.726087-02 9.990297-01 9.99307E-OI 

- 1 . 776827-01 4. 166197 + 01 3.69140E-02 9.99049E-01 9.993217-01 

-2.868/67-01 4.149767+01 3.G7/807-02 9.990567-01 9.993267-01 

-3.12191C-01 4.13183E+01 3.661917-02 9.990637-01 9. 993307-01 

-3 4.113 13^+01 3.64696E-02 9.99073F-01 9.993377~01 

-4.130S9F-01 4.093127+01 3.62769E-02 9.99083E-01 9.99346L-01 

-1.270137-01 4.078027+01 3.614216-02 9.99059E-01 9.993496-01 

-2.376 337-0 1 4. I 78306 + 01 3. 703646-02 9.990436-01 9.993167-01 

-3.98662F-02 4.138477+01 3.66779E-02 9.990596-01 9.993286-01 

-1.135067-01 4.122036+01 3.65322E-02 9.990696-01 9.993356-01 

-1.704276-01 4.106106+01 3.639096-02 9.990766-01 9.99340C-0I 

-2.130846-01 4.094197+01 3.620546-0? 9.99079F-0I 9.993427-01 

-3-1^50807-01 4.002807 + 01 3.610527-02 9. 990066-01 9.993476-01 

-1.427127-01 4.071307+01 3.608256-02 9.990697-01 9. 993496-01 

-7. 5t>7f)8L-01 3.926107+01 3.47953F-02 9.99I56F-01 9.99397F-01 

-a. 375127-01 3.39757F+01 3.454246-02 9.991667-01 9.994047-01 

-4.709706-01 3.R4655E+01 3.409016-02 9.991696-01 9.994216-01 

-5.571597-01 3.B2782E+01 3.392407-02 9.991966-01 9.99426F-01 

-3.283137-01 3.838117+01 3.401716-02 9.991936-01 9.994236-01 

-1.3cl315E-0l 3. 660797+01 3.42072F-02 9.99179F-QI 9.994147-01 

-6.220706-01 3.911677+01 3-466746-02 9.991596-01 9.903997-01 

-7.379295-01 4.040007+01 3.500506-02 9.991066-01 9.993617-01 

-l.ri766E+00 4.01548F+01 3.576506-02 9.991066-01 9.99361F-0I 

-0.19441F-OI 3.997506+01 3.542996-02 9.991237-01 9.993736-01 

-6.862897-01 1.978346+01 3.525846-02 9.991336-01 9.993006-01 

-4. 412497-01 3.9694GE + 0I 3.517966,-02 9.991366-01 9.993016-01 

-4. 644106-01 3.975436+01 3.523266-02 9.991336-01 9.993807-01 

1. 614716-01 3.977J0C+0I 3.524926-02 9.991336-01 9.993006-01 

-1.006047+00 4.205306+01 3.727116-02 9.990296-01 9.993076-01 

-7.36780F.-Ol 4.179416 + 01 3.704096-02 9.990436-01 9.993167-01 

-9.6650St-Ol 4.166646+01 3.692776-02 9.990497-01 9.993217-01 

-9. 450317-01 4.146137 + 01 3.674596-02 9.990566-01 9.993266-01 

-9.974246-01 4.121637+01 3.652876-02 9.990697-01 9.993356-01 

-9 . 798,>3F-0 1 4. 096356 + 01 3.630456-02 9.990796-01 9.993426-01 

-5.179667-01 4.057566+01 3.596076-02 9.990966-01 9.993547-01 

-i. 539086+00 4.384347+01 3.805777-0? 9.909466-01 9.992477-01 

-0.303376-01 4.3352^6+01 3.842286-02 9.98969F-01 9.992647-01 

-I f I 1 v5537 + 00 4.339067+01 3. 045626-02 9.989666-01 9.992617-01 

-1.133077+00 4.312797+01 3.022336-02 9.93979F-01 9.992717-01 

-1.408127+00 4.272246+01 3.786307-02 0.989996-01 9*992057-01 

-1.496047+00 4.212536+01 3.733456-02 9.990266-01 9.993047-01 

-0. 850837-01 4.14005F+01 3.669196-02 9.99059E-0I 9.9932OF-0I 

-1.424307+00 4.42<U3C+01 3.923477-02 9.989236-01 9.992306-01 

-9.O2710L-01 4.375716+01 3.070116-02 9.909496-01 9.992497-01 

-9.5059Or-ol 4.36995F+01 3.073OIF-O2 9.989536-01 9.992526-01 

-1.052007+00 4.33958F+01 3.046086-02 9.909666-01 9.992616-01 

-1. 243096+00 4.295817+01 3.007206-02 9.989096-01 9.992706-01 

-1.615297+00 4. 235357+01 3.754126-02 9.990167-01 9.992976-01 

-9.01247F-01 4.167656+01 3.693666-02 9.990466-01 9.993107-01 

1.54960F-OI 3.603846+01 3.19386E-02 9.992896-01 9.994926-01 


1 




FIG. 1. SCROLL SCHEMATIC SHOWING THE COORDINATES 
AND THE SPLITTING SURFACES. 
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FIG. 2. STRUCTURE OF THE PROGRAM. 
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FIG. 3. MODAL PLACEMENT ON THE SCROLL SURFACE ABC. 
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FIG, 6, CONTOURS OF COMPUTED CIRCUMFERENTIAL VELOCITY COMPONENT. 





APPENDIX A 


DERIVATION OF THE ELEMENT EQUATIONS 


The Simple four node linear tetrahedral element is used 
in the present analysis. In this case the unknown velocity 
potential field in an element is given by the following 
equation: 

e i 

<{> = I <i>. N. (A.l) 

i=l ^ ^ 

where (f)^ are the unknown velocity potentials at the elements 
four nodes , and 

N^ are the linear interpolation functions (i=l,2,3,4). 

The element equation is given by equation (5) which is 
rewritten below: 

/ VN^*(pV[N] {(})}®)dV = / p 1^ dA (A. 2) 

V® A® 

The above equation will be written in a different form as 
follows : 

[K]® {(j)}® = {P}® (A. 3) 

The Stiffness Matrix 

The coefficients of the element stiffness matrix can be 

shown [5, 7] to be given by: 


3N. 3N. 3N 3N. 3N . 3N . 

^i j ~ ^ 3x 3x^ 3y 3y^ ^ 3z 9z^ ^ ^ 


(A. 4) 


An element matrix [B]^z is defined to contain the geometric 
properties of the element as follows : 
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1 Yl 

1 Xj yj Zj 

1 Xj yj Z3 

1 y^ z^ 


(A. 5) 


In the following derivations, it will be assimed that the 
local numbering of the volume element nodes is such that the 
first three are numbered counterclockwise in ascending order 
when viewed from the fourth node [7] . 

Using linear interpolation functions, 


N 


^ = a^x + h^Y + c^2 + d^ 


(A. 6) 


it can be shown that the derivatives of the interpolation 
functions are given by [7] , 


8N 

S3T = ^i = ^ cofactor (B^_ ) 


3N 

■r — = b = — — cofactor (B ) 

9y 1 


3N. , 

•77 — - = c. = cofactor (B ) 

3z 1 _„e z. 

6V 1 


(A. 7) 


The element stiffness matrix is then given by 


a. a. +b.b. +c.c. 
K., = p(-^J LJ. 


1 : 


36 V" 


(A. 8) 


where 


V® = -| det [B]® 


(A. 9) 
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Equation (A. 8) is used to build up the global matrix of coeffi- 
cients in subroutine MATVEC. 

The load vector surface integrals on the right hand side 
of equation (A. 2) are only evaluated for the surfaces of the 
element that lie on the external boundaries with mass flow 
across. The contribution of the surface element to the load 
vector corresponding to each of its 3 nodes is given by: 



where f is the mass flux into the element surface whose area 
IS equal to A®. Equation (A. 10) is used to build up the 
components of the load vector in subroutine MATVEC. 


The Flow Velocity 

The velocity in the linear tetrahedral elements 
can be expressed in terms of the derivative of the inter- 
plation functions and velocity potential at element nodes 
as follows: 


4 9N. 


V 


X 


= ^ ^i = ^i'f’i 


1=1 


4 9N 

V = y ^ (j) . = b (f) . 


4 3N. 


V 


= y -5 — - 4) . = c (j) . 

3z 1 11 


i=l 


(A. 11) 


where the repeated indices imply summation over i=l to 4. 
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APPENDIX B 


THE MATCHING CONDITIONS AT THE SPLITTING SURFACES 


The splitting surfaces are introduced in the solution 
domain in order to be able to handle the mathematically multi- 
valued velocity potential in the numerical solution. They can 
be arbitrarily placed in the flow field extending from the 
point where,two streams of different histories join at the scroll 
lip^ to the exit station. Nodes with different global numbers, 
but which have the same coordinates are placed opposite each 
other on the two splitting surfaces. The velocity potentials 
at each two nodes will generally be different, thus allowing 
for a velocity potential discontinuity in the numerical 
solution. Conditions are introduced in the equations to match 
the mass flow across the splitting boundaries and the difference 
between the velocity potential values of the opposite nodes 
on the splitting surfaces. 

The first condition imposed in the numerical solution 
is the matching of the weighted mass flow across the splitting 
surfaces. Across an element's surface, s®, on the splitting 
boundary, the weighted mass flow is given by: 




d A^ 


(B.l) 


where 3(|)®/3n is derivative of the velocity potential normal 

to the surface area A^. This derivative can be expressed in 

terms of the x, y and z derivatives of (j) and the direction 

cosines a, 3 and y of the normal to the surface A as follows: 

s 


d<t>^ 

9n 


4 3N. 4 ^ 


3N. 4 3N. 


3y 


(B.2) 


Substituting equation (B.2) into (B.l) and using equation (A. 7) 
one obtains the following expression for the contribution of the 
weighted mass flow across the splitting surface to each of the 
opposite surface nodes. 
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E 


1 

3 


4 

I p(aa^ + 6b^ + YC^) 


4>.A® 

s 


(B.3) 


The direction cosines a, S and y of the normal to the surface of 
the element on the splitting boundary can be expressed in terms 
of the coordinates of its three nodal points. If the first 
three nodes of a volume element lie on the splitting surface 
and the fourth is an internal node, it can be shown that the 
direction cosines a, 3 and y proportional to cofactor , 


cofactor By^ and cofactor b^^ respectively, 
constant is determined from 


The proportionality 


2 2 2 

a + B + y = 1 . (B.4) 

The second condition is imposed to match the velocity potential 
difference between each two opposite nodes on the splitting 
surfaces. In the following derivations and discussions, 
it will be assumed that if the nodes of the volime element on 
one of the splitting boundaries are i, j, k, etc., the opposite 
nodes on the other splitting boundary are i', 3 ', k', etc. The 
mass flow matching condition is introduced by adding nig/ 3 to 
the global matrix coefficients in the columns corresponding to 
the nodes (i = 1,2, 3, 4) in the rows (i = l',2',3',4') in sub- 
routine SPLIT. The same process is repeated reversing the roles 
of the primed and unprimed sides. 

The second matching conditions insures the same difference 
between the velocity potential values of the opposite nodes 
on the splitting surface. 

I t 

• (3.5) 

This condition is introduced by modifying the coefficients of the 
row corresponding to the node i on one side of the splitting 
boundaries. This is accomplished by replacing the diagonal 
coefficient and the coefficient in column j ' by a very large 
number and those in columns i' and j by equal but negative 
large numbers. 
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APPENDIX C 


NOMENCLATURE 

A Surface area. 

B The element matrix. 

K Matrix of coefficients or stiffness matrix, 

m mass flow rate. 

N Interpolation function. 

P Load vector 

R Gas constant. 

T Temperature. 

x,y,z Cartesian coordinates 

a,6,Y Direction cosines of the normal to a surface 

Y Specific heat ratio, 

p Gas density. 

(J) Velocity potential 

Superscripts 

e Refers to the volume element. 

Subscripts 

i,j Refers to the nodes of an element or the indices 

of the corresponding stiffness matrix, 
s Refers to the splitting surfaces. 
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End of Document 



